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PREFACE 


The  objective  of  this  research  project  is  to  understand  phenomena  of  seismic  wave  generation 
and  propagation  that  affect  regional  and  teleseismic  seismograms  used  for  nuclear  test  mon¬ 
itoring.  Two  phenomena  were  investigated  as  part  of  this  project:  scattering  by  irregular 
interfaces,  and  the  effect  of  non-sphericity  of  a  cavity  containing  an  explosion  on  its  seismic 
wave  radiation.  Preliminary  results  on  the  latter  topic  have  been  reported  in  last  year’s  an¬ 
nual  techincal  report  and  at  the  15th  Annual  Seismic  Research  Symposium  at  Thornwood, 
New  York.  Our  work  on  modeling  the  effects  of  non-spherical  tunnels  and  cavities  continues 
under  an  Air  Force  Office  of  Scientific  Rersearch  grant  in  progress,  and  final  results  of  this 
work  will  be  reported  in  future  reports  and  publications. 

This  final  report  contains  results  of  our  invetigations  of  interface  scattering.  Under  an 
earlier  contract  we  reported  the  development  of  theoretical  methods  for  modeling  scattering 
from  highly  irregular  interfaces.  Here  we  report  on  the  application  of  these  methods  to 
laboratory  and  field  data.  The  first  section  of  the  report  is  a  preprint  of  a  paper  submitted 
to  the  Journal  of  the  Acoustical  Society  of  America.  It  describes  the  results  of  ultrasonic 
laboratory  experiments  involving  backscattering  of  acoustic  waves  from  a  glass  surface  etched 
to  produce  a  highly  irregular  3-D  interface.  The  laboratory  measurements  are  compared  to 
numerical  simulations  performed  with  boundary  element  modeling.  The  second  part  of  this 
report  is  a  study  of  P  wave  coda  observed  at  small  aperture  arrays  in  Scandinavia  and  New 
England  from  events  at  regional  distances.  The  frequency-wavenumber  characteristics  of 
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the  coda  are  analyzed  and  interpreted  in  terms  of  scattering  from  an  irregular  crust-mantle 
boundary  or  interfaces  within  the  crust.  This  part  of  the  report  has  been  submitted  for 
publication  in  the  Bulletin  of  the  Seismological  Society  of  America. 
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EXPERIMENTAL  STUDY  OF  SCATTERING  FROM  A 


HIGHLY  IRREGULAR,  ACOUSTIC-ELASTIC  INTERFACE 
Summary 

In  this  study,  we  experimentally  and  numerically  develop  statistical  models  for  the  scattering 
of  an  acoustic  P  wave  which  is  incident  on  a  highly  irregular,  random  acoustic-elastic  interface 
to  determine  the  general  nature  of  reflected  energy.  We  then  elucidate  whether  or  not 
enhanced  backscattering,  already  identified  numerically  for  acoustic  (SH)  and  fully  elastic 
media  (P-SV),  occurs.  Numerically,  the  problem  is  solved  in  two  dimensions  by  coupling  the 
representation  theorem  for  an  elastic  medium  with  the  extinction  theorem.  Exact  integral 
expressions  for  the  scattered  pressure  in  the  acoustic  medium  are  then  obtained,  which 
include  all  converted  and  all  multiply  scattered  waves  at  the  boundary.  Experimentally,  a 
glass  etching  process  using  photoresist  templates  with  Gaussian  statistics  allowed  for  the 
generation  of  characterized  interface  irregularities.  Experiments  were  performed  on  a  glass 
surface  with  an  rms  slope  of  30  degrees,  for  the  case  of  an  incident  wavelength  with  a  size  on 
the  same  order  as  the  interface  irregularities.  The  numerical  models  predict  an  enhancement 
of  energy  difiracted  back  towards  the  source,  and  results  obtained  in  our  in-house  ultrasonic 
laboratory,  strongly  support  the  presence  of  this  retroreflective  energy.  In  terms  of  general 
scattering,  we  find  that,  at  smaller  incident  angles  (relative  to  vertical),  the  2-D  numerical 
results  can  give  insight  into  the  3-D  experimentally  observed  scattering  over  most  scattering 
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angles.  However,  at  larger  incident  angles,  fundamental  differences  between  2-D  and  3-D 
scattering  may  exist. 

Introduction 

In  laboratory  experiments,  a  lack  of  control  over  the  statistical  parameters  of  a  given  ran¬ 
dom  model  can  easily  produce  ambiguous  results.  In  the  case  of  irregular  interfaces,  the 
height  probability  distribution  and  the  correlation  lengths  of  the  interface  may  be  poorly 
constrained,  the  interface  statistics  may  show  nonstationarity,  and  the  interface  may  contain 
a  wide  variety  of  length  scales.  Each  of  these  experimental  uncertainties  makes  comparisons 
with  numerical  models  difficult,  if  not  impossible.  It  is  the  goal  of  this  study  to  physically 
fabricate  a  random  interface  which  is  stationary  in  space,  with  both  a  simple  probability 
distribution  in  height  and  a  simple  transverse  correlation  function.  Experimental  results  can 
then  be  easily  compared  with  the  corresponding  numerical  results. 

The  accurate  physical  generation  of  a  Gaussian  surface  is  important  experimentally,  since 
Gaussian  interfaces  are  mathematically  convenient  and  have  been  widely  used  up  to  this  point 
in  scattering  studies.  Many  theoretical  formulations  in  the  literature  apply  the  simple  prop)- 
erties  of  a  Gaussian  correlation  function  to  heterogeneities.  Examples  can  be  found  in  Prange 
and  Toksoz  (1990),  Knopoff  and  Hudson  (1964,  1967),  Haddon  (1978),  and  Kuperman  and 
Schmidt  (1989).  However,  a  Gaussian  autocorrelation  function  is  a  somewhat  unrealistic 
model  for  many  regions  of  the  earth  since  a  Gaussian  function  is  continuously  differentiable 
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and  has  a  rapid  decay  in  spectral  amplitude.  This  continuous  differentiability  and  narrow 
wavenumber  spectrum  result  in  a  very  smooth  model  with  one  dominant  length  scale.  More 
general  fractal  processes,  which  exhibit  fluctuations  on  all  length  scales,  are  likely  a  more 
realistic  model  of  true  Earth  structure.  The  most  commonly  utilized  fractal  representation 
is  an  exponential  covariance  function.  In  this  case  the  model  is  continuous  but  not  differ¬ 
entiable,  making  it  a  far  rougher  model  than  that  given  by  the  Gaussian  function.  The 
exponential  covariance  function  has  been  utilized  extensively,  and  in  many  instances  gives  a 
good  description  of  physical  properties  (e.g.,  Wu  and  Aki,  1985;  Frankel  and  Clayton,  1986). 
Recently,  Goff  and  Jordan  (1988)  have  generalized  the  exponential  function.  In  this  case, 
seafloor  topography  was  modeled  as  a  two-point  covariance  function  with  five  free  parame¬ 
ters  which  describe  the  amplitude,  orientation,  characteristic  wavenumbers,  and  Hausdorff 
(fractal)  dimension  of  the  topography.  These  self-affine  surfaces,  of  which  the  exponential 
covariance  function  is  a  special  case,  have  been  shown  to  give  good  first  order  stochastic 
descriptions  of  seafloor  morphology. 

As  a  first  step  at  modeling  scattering  by  irregular  interfaces,  we  fabricate  and  physically 
model  an  irregular  interface  with  a  Gaussian  correlation  function.  Upon  evaluating  the 
effectiveness  of  this  experimental  approach  for  the  smoother  Gaussian  surface,  more  general, 
self-affine  interface  models  may  be  proposed.  In  this  study,  the  statistical  parameters  of 
the  interface  are  chosen  so  that  the  incident  wavelength  has  the  same  length  scale  as  the 
correlation  length  of  the  irregularities.  In  addition,  the  rms  slope  of  the  interface  is  chosen  to 
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be  large,  placing  the  model  in  a  regime  where  approximate  techniques  such  as  the  Kirchhoff, 
Born,  and  the  geometrical  ray  approaches  break  down  and  multiple  scattering  mechanisms 
such  as  “enhanced  backscattering”  and  “shadowing”  play  dominant  roles. 

What  is  “enhanced  backscattering”  from  an  acoustic-elastic  halfspace?  By  definition, 
enhanced  backscattering  or  “retroreflectance”  is  the  enhancement  of  energy  scattered  back 
in  the  direction  of  the  source.  In  optical  theory,  O’Donnell  and  Mendez  (1987)  were  the  first 
to  propose  the  hypothesis  that  time-reversed  paths  are  responsible  for  enhanced  backscatter¬ 
ing.  This  hypothesis  was  further  strengthened  by  Maradudin  et  al.  (1990)  who  showed  that 
retroreflectance  exists  first  for  energy  double  scattered  from  an  interface.  Further  support 
for  this  hypothesis  came  from  the  work  of  Schultz  and  Toksoz  (1993,1994),  which  showed 
numerically  that  retroreflectance,  which  exists  in  the  general  case  of  seismic  scattering,  is 
consistent  with  the  time-reversed  path  hypothesis.  P-SV  scattering  gave  the  strongest  sup¬ 
port,  as  the  hypothesis  could  be  tested  on  various  wave  conversions  at  the  interface.  As 
predicted,  energy  enhancement  was  clearly  observed  on  P-to-P  and  S-to-S  scattering  and 
not  on  P-to-S  and  S-to-P  scattering. 

The  idea  of  time-reversed  paths  is  easily  extended  to  the  acoustic-elastic  case.  Take  for 
instance  the  peak-valley  sequence  shown  in  Figure  1.  If  an  incident  P-wave,  shown  by  the 
solid  line,  diffracts  from  point  1  it  will  propagate  as  a  P-wave  to  point  2  and  then  diffract 
at  some  angle  into  the  upper  medium  again  as  a  P-wave.  For  most  waves  traveling  away 
from  the  interface,  enhancement  will  not  occur.  However,  if  the  diffracted  wave  travels 
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directly  back  towards  the  source,  a  plane  wave  component  of  incident  P-wave  can  be  found 
traveling  exactly  in  the  reverse  direction:  propagating  from  point  2  to  point  1,  and  traveling 
back  towards  the  source  as  shown  by  the  dashed  line.  In  this  case,  the  time-reversed  path 
interferes  constructively  with  the  forward  path  and  contributes  additional  energy  towards 
the  source  resulting  in  enhanced  backscattering.  Using  a  simple  phase  argument  based  on 
each  of  the  studies  cited  above,  some  properties  of  enhanced  backscattering  can  be  derived. 
The  peak  width  can  be  written  as  A6s  =  j,  where  AOg  is  the  angular  width  of  the  peak,  A 
is  the  incident  wavelength,  and  I  is  the  mean  free  path  of  the  interface  or,  in  other  words, 
the  average  distance  a  wave  propagates  between  points  1  and  2  along  the  interface. 

Adapting  this  phase  approach,  other  path  geometries  may  also  contribute  to  enhanced 
backscattering.  For  example,  if  a  wave  encountering  a  3-D  irregular  interface  multiply  scat¬ 
ters  from  several  points  outside  the  vertical  source-receiver  plane  and  then  sends  energy 
back  in  the  direction  of  the  source,  then  a  time-reversed  path  may  be  found  which  also 
sends  energy  back  towards  the  source.  In  the  same  manner,  many  multiply  scattered  paths, 
sending  energy  back  towards  the  source,  can  be  found  in  the  acoustic-elastic  case.  However, 
as  a  result  of  energy  loss  with  each  diffraction  from  the  interface,  due  to  both  transmission 
through  the  interface  and  additional  spreading,  it  seems  reasonable  that  the  double-scattered 
paths,  both  in  and  out  of  the  source-receiver  (incident)  plane,  will  contribute  the  majority 
of  retrorefiective  energy. 

Experimentally,  there  has  been  little  investigation  into  coherent  backscattering  in  the 
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acoustic  case.  However,  for  disordered  acoustical  systems,  the  recent  study  of  Bayer  and 
Niederdrank  (1993)  has  shown  that  strongly  heterogeneous  media,  such  as  random  distribu¬ 
tions  of  brass  rods  (2-D)  and  random  distributions  of  gravel  (3-D)  can  create  weak  localization 
effects.  Coherent  backscattering  was  studied  in  a  narrow  region  about  the  backscattering 
direction  where  a  simple  diffusion  model  predicted  the  results  well.  Theoretical  work  has 
been  done  on  the  localization  of  acoustic  waves,  where  Kirkpatrick  (1985)  has  shown  that 
all  states  are  localized  for  dimension,  d  <2.  Although,  this  work  may  have  implications  for 
seismic  scattering  in  an  inhomogeneous  crust,  Condat  and  Kirkpatrick  (1987)  have  shown 
that  the  locahzation  of  waves  in  an  acoustic  medium  is  very  difficult  to  achieve  unless  ex¬ 
tremely  strong  scatterers  are  present.  Localization  resulting  from  crustal  heterogeneities, 
therefore,  is  likely  far  smaller  than  that  resulting  from  interfaces. 

This  study  is  organized  as  follows.  The  first  section  briefly  summarizes  the  numerical  for¬ 
mulation  used  to  model  scattering  from  a  randomly  irregular  acoustic-elastic  interface.  This 
Somigliana  identity  approach  is  an  extension  of  the  work  of  Schultz  and  Toksoz  (1993,1994). 
In  the  second  section,  the  physical  construction  of  the  characterized  glass  surface  is  discussed 
and  the  ultrasonic  apparatus  for  measuring  the  amplitude  distribution  of  energy  scattered  by 
the  interface  is  described.  The  third  section  directly  compares  3-D  experimental  results  ob¬ 
tained  in  our  in-house  ultrasonic  water  tank  with  the  corresponding  2-D  numerical  results. 
Amplitude  distributions  are  presented  in  detail  and  scattering  mechanisms  are  proposed. 
Finally,  we  investigate  the  differences  between  the  2-D  synthetic  and  the  3-D  experimental 
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data  and  discuss  the  implications  for  both  2-D  and  3-D  scattering  mechanisms. 

Theory 

The  numerical  approach  and  notation  in  this  section  follows  that  of  Schultz  and  Toksoz  (1993, 
1994).  Since  the  approach  here  is  very  similar  to  the  SH  and  P-SV  cases,  we  give  only  a 
brief  outline  of  the  theoretical  approach.  We  first  express  the  total  scattered  displacement  at 
any  point  within  two  volumes  of  elastic  material  with  the  Somigliana  representation  theorem 
(e.g.,  Aki  and  Richards  1980).  Simplifying  this  theorem  to  a  2-D  case  gives  the  following 
integral  equations 

H{l]u^^\x)  =  ^/^')(2)Gg(x;2)dF(77) 

+  [dS(z'){[c\%ghj(x')dG^^l{x;^  (1) 

J  s 

-Gm(x;x')rm(u<')(x'),a)}, 

where  gradients  are  zero  in  the  0:2 -direction.  Following  the  notation  of  Schultz  and  Toksoz 
(1993),  r«(x)  is  the  traction  vector  along  the  interface  in  both  the  fluid  {I  =  /)  and  the 
solid  {I  =  s),  and  we  have  assumed  all  surfaces  to  be  far  enough  away  so  that  only  the 
surface,  S'(x),  separating  the  two  volumes,  contributes  to  the  final  displacement.  Gnp{x;x') 
is  a  Green’s  function  giving  the  n-component  of  displacement  at  x  resulting  from  a  point 
force  in  the  p-direction  at  x  ,  Cjjpg  is  the  elasticity  tensor,  and  H[i]  is  a  function  that  takes  a 
value  of  0  or  1  depending  on  whether  the  point  x  lies  outside  or  inside  the  volume  of  interest. 
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I,  respectively.  We  assume  that  the  media  are  homogeneous  and  isotropic. 


In  this  work  the  upper  medium  is  acoustic,  supporting  propagation  of  only  dilatational 
waves,  while  the  lower  medium  is  taken  to  be  elastic.  The  boundary  conditions  for  the 
resulting  acoustic-elastic  boundary,  can  be  written  in  the  general  form 

n-u(^^(x)U3=^(^,)  =  n-u('’)(x)U=C(^i), 

T(^Hx)|.3=C(xO  =  T^^Hx)|x3=c(xO>  (2) 

n  X  T(")(x)|a,3=:C(xx)  =  0, 

where  the  surface  profile  function  is  taken  to  be  0:3  =  (^(xi)  and  the  unit  normal  vector  along 
the  interface  can  be  written  as 

n  =  (-C'(x,),l)[l  +  {C'(a:i)n-i  (3) 

The  first  boundary  condition  in  (2)  represents  the  continuity  of  normal  displacement  and 
the  other  two  conditions  together  represent  the  continuity  of  normal  stress. 

Taking  our  volume  of  interest  to  be  the  upper  acoustic  medium,  placing  the  incident 
wave  in  the  acoustic  medium,  substituting  the  final  form  of  the  boundary  conditions  (2), 
and  letting  X3  — >  ^'*’(3:1),  the  final  set  of  coupled  integral  equations  can  be  written  as 

S{x)  =  S^^\x)incid  (4) 

/-f  00 

da:'i[5(x')T(^nx|x')  -  Di-^)(x|x')T>n(x')], 

-00 
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in  the  acoustic  medium,  and 


0  =  -  T"  dx[lUi{^)T^%\x')  -  ^Di’\x\^)S{^)], 

J—OO  H'  ' 


(5) 


in  the  elastic  medium,  where  we  have  defined 


Di^\x\x')  =  IJ-^"''G^nii^\x')npUs=cixx)- 


(6) 


in  the  respective  media.  Now  the  unknown  source  strength  functions,  which  we  eventually 
solve  for,  can  be  expressed  as  a  function  of  Xi  alone 


S(l,)  =  (7) 


where  6{x)  =  Uk,k{'^)  is  the  dilatational  parameter,  A  and  (i  are  the  Lame  parameters  for  the 
medium,  U(x)  is  the  displacement  along  the  interface,  where  should  not  be  confused 

with  the  surface  function  referred  to  previously.  The  pressure  term  has  been  normalized 
with  respect  to  A^-^^  to  ensure  that  the  final  numerical  system  is  well  conditioned. 

We  now  approximate  the  transducer  source  as  a  narrow  Gaussian  beam  source  (Schultz 
and  Toksoz,  1993).  This  beam  source  acts  to  dramatically  reduces  the  computational  demand 
of  this  numerical  approach,  since  only  a  small  portion  of  the  interface  is  excited  by  the 
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incident  beam  and  the  length  of  integration  along  the  acoustic  boundary  is  reduced.  The 
half  width  of  the  incident  beam  is  w  =  h  cos  do,  where  h  is  the  half-width  of  the  incident 
beam  as  projected  on  the  plane  0:3  =  0.  The  scattered  field  in  the  acoustic  medium  can 
now  be  expressed  completely  in  terms  of  the  unknown  source  functions.  Substituting  the 
Cartesian  form  of  the  Green’s  functions  for  a  homogeneous  halfspace  (Bouchon  and  Aki, 
1977),  the  scattered  pressure  at  any  point,  X3  >  C{^i)max,  iri  the  fluid  can  be  decomposed 
into  a  summation  of  plane  waves. 

An  approximation  to  the  Fourier  reflection  coefl&cient  can  now  be  written  in  terms  of 
an  amplitude  coefficient  by  normalizing  the  scattered  wave  amplitude  coefficient  by  the 
amplitude  of  the  incident  plane  wave  component  having  a  wavenumber  vector  corresponding 
to  the  angle  of  incidence.  The  reflection  coefficient  can  then  be  expressed  as 


where 


R{kuj)  = 


2y/Tik\^^'w' 


(8) 


fpi^s)  =  [  da;i[i5(a:i)A:p^(sin0,C(2^i)  “  COS0,)  (9) 

J—OO 

-ki^^^{C{x[}D,{x[)  -  £)3(a:;))]e-fel''(sin0,xi+cos^,C(xi))^ 


which  is  comparable  in  amplitude  to  the  Fourier  reflection  coefl&cient  calculated  for  a  single 
incident  plane  wave.  This  normalization  is  diflPerent  from  the  Diflferential  Reflection  Coeffi¬ 
cient  implemented  by  Schultz  and  Toksoz  (1993,1994).  Note  that  we  let  k  =  fc^^sin^s  and 
kz  =  kt^^  cos  6  s- 
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These  integral  equations  can  now  be  solved  numerically.  The  solution  over  an  acoustic- 
elastic  interface  can  be  expressed  completely  as  a  combination  of  the  Green’s  functions  for 
the  P-SV  and  SH  cases  (Schultz  and  Toksoz,  1993,  1994),  where  the  shear  velocity  of  the  SH 
wave  is  changed  to  the  P-wave  velocity  of  the  acoustic  medium,  so  as  to  reflect  the  acoustic 
Green’s  function.  The  flnal  coupled  integral  equations  are  then  transformed  to  a  coupled  set 
of  matrix  equations  and  solved  using  LU  decomposition  or  biconjugate  gradient  techniques. 

In  this  paper  two  types  of  interfaces  are  modeled  numerically.  Both  interfaces  have  a 
Gaussian  distribution  about  the  mean,  where  6^  =  (C^(a:i))  is  the  mean-square  departure  of 
the  surface  from  flatness.  The  first  interface  studied  has  a  correlation  function 

W{\xi  -  xil)  =  r2(C(xi)C(xi)),  (10) 

described  by  the  Gaussian  function,  lP’(|a;i|)  =  exp{—x\/ a?).  The  correlation  length,  a, 
for  a  Gaussian  interface  is  approximately  equal  to  the  average  distance  between  adjacent 
peaks  and  valleys  along  the  interface.  The  interface  can  also  be  described  in  terms  of  the 
rms  slope  of  the  surface,  0,  which  we  will  refer  to  often.  This  rms  slope  can  be  written 
as  0  =  tan~^{^).  The  second  surface  studied  has  an  exponential  correlation  function, 
kP(|a;i|)  =  exp{-xi/a). 

Averaging  over  an  ensemble  of  realizations  of  these  surfaces,  we  display  the  final  scattered 
pressure  as  a  statistical  average  that  follows  the  approach  of  Schultz  and  Toksoz  (1993,1994). 
The  total  mean  squared  contribution  to  the  reflection  coefficient  (RC)  can  then  be  written 
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as 


This  gives  the  average  square  pressure  reflected  into  the  upper  medium  as  a  function  of  the 
scattering  angle,  6s,  given  one  incident  beam  angle,  ^o-  The  square  root  of  the  RC  is  used  for 
comparison  with  experimentally  recorded  amplitudes.  Given  this  representation,  a  Gaussian 
beam  which  is  perfectly  reflected  by  a  plane  interface  gives  a  maximum  reflection  coefficient 
value  of  one. 


Experimental  Procedure 

The  experimental  approach  consists  of  submerging  a  solid  elastic  model,  in  this  case  a  glass 
block,  into  an  experimental  water  tank,  thus  creating  an  acoustic-elastic  interface.  More 
precisely,  this  experiment  involved  the  generation  of  an  interface  with  a  characterized  random 
surface  and  the  design  of  a  measurement  configuration  that  could  accurately,  to  within  a 
fraction  of  a  degree,  measure  various  realizations  of  the  fluid-glass  boundary. 

Random  Interface  Generation 

Fabricating  the  randomly  irregular  surface  was  the  most  challenging  part  of  this  project. 
Numerous  irregular  surfaces  were  designed.  Models  ranged  from  irregular  distributions  of 
glass  beads  to  roughened  cement  surfaces.  In  addition,  random  distributions  of  gravel  held 
together  by  epoxy  were  tested  along  with  naturally  rough  granite  and  sandstone  surfaces. 
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Unfortunately,  these  models  either  did  not  give  proper  control  over  statistical  parameters  or 
were  extremely  heterogeneous  at  the  ultrasonic  level,  making  comparisons  with  numerical 
models  very  difficult.  After  much  experimentation  the  most  promising  approach  became  the 
fabrication  of  a  random  glass  surface  using  an  etching  procedure. 

The  irregular  glass  surface  was  designed  using  a  solid  glass  block  and  a  standard  glass 
etching  process.  First  a  cylindrical  glass  block  with  a  diameter  of  19.5  cm  and  a  height 
of  7.5  cm  was  cast  using  a  graphite  mold.  Then  to  achieve  the  desired  random  interface, 
the  Gaussian  surface  described  in  the  previous  section  was  first  numerically  generated.  The 
transverse  correlation  length,  a,  and  the  standard  deviation  of  the  height,  5,  of  the  interface 
were  specified  as  1  mm  and  .71  mm,  respectively,  giving  an  rms  slope  of  45°.  After  its 
numerical  generation,  the  Gaussian  surface  was  discretized  into  six  individual  depth  levels, 
with  each  level’s  thickness  equal  to  one  standard  deviation  of  the  interface  height,  6.  The 
templates,  shown  in  Figure  2,  were  successively  glued  to  the  surface  and  each  template 
was  exposed  to  high  velocity  sand  particles  normally  incident  on  the  surface.  As  a  result, 
the  desired  interface  was  etched  in  six  discrete  depth  intervals.  In  general,  the  correlation 
length  of  the  surface  was  controlled  by  the  template  design  and  the  standard  deviation  of 
the  interface  was  controlled  by  the  time  that  each  template  was  exposed  to  the  sand  blast. 
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The  Scattering  Instrumentation 


Once  the  irregular  glass  surface  was  created,  an  automated  scattering  apparatus  was  used 
to  measure  the  scattering  properties  of  the  interface.  Two  different  flat-sub-bottomed  trans¬ 
ducers  were  used  to  create  a  beam  source.  They  consisted  of  a  Panametrics  12.7  mm  di¬ 
ameter  transducer  (1.5  MHz,  A  =  1.0  mm  in  water)  and  a  Panametrics  25.4  mm  diameter 
transducer  (0.5  MHz,  A  =  3.0  mm  in  water).  The  detectors,  which  were  also  Panametrics 
flat-sub-bottomed  transducers,  consisted  of  a  6.4  mm  diameter  transducer  (1.5  MHz,  A  =  1.0 
mm  in  water)  and  a  12.7  mm  diameter  transducer  (0.5  MHz,  A  =  3.0  mm  in  water),  respec¬ 
tively.  The  detectors  were  chosen  such  that  they  were  sensitive  only  to  waves  approaching 
nearly  perpendicular  to  the  bottom  surface  of  the  transducer,  limiting  the  energy  recorded 
to  waves  approaching  in  line  with  the  transducer  axis.  Given  the  source  parameters,  the 
resulting  source  radiation  pattern  was  a  beam  with  most  of  the  energy  traveling  in  the  for¬ 
ward  direction.  The  source  radiation  patterns  exhibited  some  slight  spreading  of  the  beam, 
although  further  tests  showed  that  this  did  not  significantly  affect  the  results. 

The  experimental  geometry  used  to  measure  the  scattered  amplitudes  is  shown  in  Fig¬ 
ure  3.  The  glass  block  was  located  at  the  center  of  the  experiment  and  the  source  and 
detector  were  then  stepped  in  a  semicircle  about  an  axis  of  rotation.  This  axis  extended  lat¬ 
erally  along  the  irregular  fiuid-glass  interface  and  perpendicular  to  the  source-receiver  plane. 
In  each  experiment  the  source  was  placed  at  a  constant  incident  angle,  ^o,  and  a  constant  .35 
m  distance  from  the  axis  of  rotation.  The  recording  angle  was  then  controlled  by  mounting 
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the  detector  on  a  motor-driven,  rotating  arm  that  held  the  detector  .30  m  from  the  given 
axis  of  rotation.  The  arm  was  then  rotated  in  0.9°  steps  about  this  axis  of  rotation.  There¬ 
fore,  the  recorded  energy  represents  scattering  in  the  plane  of  incidence.  Unfortunately,  the 
detector  occludes  the  source  when  it  is  near  the  backscattering  position,  resulting  in  a  loss 
of  3°  to  6°  of  measurements  near  the  source  position.  These  data  points  are  not  plotted. 

The  final  measurement  is  the  mean  scattered  pressure  as  a  function  of  scattering  angle, 
given  a  fixed  angle  of  incidence.  The  scattered  pressure  is  measured  at  specific  frequencies  by 
introducing  directly  via  the  transducer  source  a  monochromatic  wave  of  a  given  frequency. 
The  continuous  wave  is  approximated  by  a  finite  sinusoid  ranging  between  35  and  100  cycles 
and  the  final  constant  amplitude  of  the  scattered  pressure  is  recorded.  Since  for  a  single 
realization,  the  waves  scattered  from  an  irregular  interface  exhibit  large  fluctuations  in  pres¬ 
sure  as  a  function  of  the  scattering  angle,  the  experimental  data  was  averaged  to  obtain  a 
final  mean  reflection  coefficient.  To  this  end,  a  simple  scheme  was  developed  to  generate 
many  realizations  using  the  same  3-D  sample.  In  general,  each  realization  was  acquired  by 
rotating  and  shifting  the  sample  in  a  prechosen  sequence.  In  general,  the  sample  was  rotated 
by  60°  staggered  steps,  followed  by  1.25  cm  lateral  shifts  of  the  block  relative  to  the  mea¬ 
suring  devices.  Since  rotating  the  surface  with  respect  to  the  incident  beam  formed  another 
scattering  geometry,  many  different  realizations  of  the  interface  were  obtained. 

As  a  result  of  our  interest  in  the  enhancement  of  energy  traveling  directly  back  towards  the 
source,  a  receiver  was  created  that  could  retrieve  energy  in  the  occluded  zone  near  the  source. 
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This  was  achieved  by  constructing  a  piezofilm  receiver,  which  had  an  impedance  similar  to 
that  of  water,  therefore,  allowing  the  source  beam  to  propagate  through  it  with  minimal 
distortion.  More  precisely,  a  four-layered,  in-parallel  piezofilm  receiver  was  constructed  to 
allow  for  a  significant  increase  in  the  signal  to  noise  ratio. 

The  data  presented  in  the  next  section  represents  the  mean  diffusely  scattered  signal,  as 
a  function  of  angle,  for  a  fixed  solid  angle  of  data  acquisition.  No  artificial  angular  factors 
are  introduced  to  the  data.  For  comparison,  the  data  is  normalized  by  the  numerical  RC 
calculated  at  normal  incidence. 


Surface  Scattering  Measurements 

In  this  section  we  present  the  average  reflection  coefficients  obtained  over  the  roughest  glass 
surface  fabricated,  thus  allowing  a  study  of  the  strongest  form  of  scattering  at  an  irregular 
interface.  Figure  4a  shows  the  target  surface  height  distribution  for  this  interface,  indepen¬ 
dent  of  lateral  position,  and  Figure  4b  shows  the  histogram  of  the  surface  height,  based  on 
surface  profilometer  measurements  of  the  actual  surface.  The  surface  measurements  have  a 
lateral  resolution  of  about  5  which  is  within  the  range  of  accuracy  required  for  resolv¬ 
ing  the  interface  statistics.  These  histograms  show  that  the  surface  follows  approximately 
a  Gaussian  probability  distribution  with  a  standard  deviation  of  about  0.6  mm,  which  is 
close  to  the  target  value  of  0.71  mm.  Figure  5a  gives  the  target  Gaussian  autocorrelation 
function,  and  Figure  5b  shows  the  the  actual  autocorrelation  function  calculated  from  pro- 
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filometer  measurements.  The  correlation  length  of  1.4  mm  is  greater  than  the  target  value 
of  1  mm.  Also  plotted  are  the  Gaussian  and  exponential  autocorrelation  functions  with  the 
same  correlation  length  as  the  data.  The  measured  autocorrelation  function  is  very  close 
to  a  Gaussian  correlation  function  at  the  more  important,  smaller  lag  distances.  At  larger 
lag  distances,  the  surface  lies  directly  between  a  Gaussian  and  an  exponential  correlation 
function. 

Figure  6a  gives  a  grayscale  plot  of  the  experimental  surface  based  on  surface  profilometer 
measurements.  Figure  6b  plots  the  surface  height  for  a  profile  taken  across  the  surface,  while 
Figure  6c  shows  a  numerically  generated  Gaussian  and  exponential  surface  given  the  same 
standard  deviation  and  correlation  length.  The  etching  process,  being  a  natural  process, 
not  surprisingly  tends  towards  a  rougher  interface  with  a  broader  power  spectrum  than  the 
targeted  Gaussian  spectrum.  However,  even  with  these  deviations,  it  is  clear  upon  referring 
to  Figure  6c  that  the  Gaussian  surface  matches  the  experimental  interface  well,  in  both  the 
observed  slopes  and  the  lateral  scale  of  the  irregularities.  Although  the  Gaussian  interface 
gives  a  good  fit  to  the  experimental  data,  out  of  interest  we  shall  also  include  the  numerical 
results  for  an  exponential  surface. 

Based  on  the  measurements  above,  the  final  glass  interface  has  approximately  a  30°  rms 
slope.  The  slopes  of  this  interface  are  quite  steep  and  the  impedance  contrast  at  the  fluid- 
glass  interface  is  large  as  the  glass  interface  has  properties  very  similar  to  those  of  an  igneous 
material.  As  a  result,  multiple-scattering  and  shadowing  effects  can  play  a  significant  role 


17 


at  both  small  and  large  incident  angles,  and  approximate  linear  theories,  such  as  the  Kirch- 
hoff  and  Born  approaches,  break  down.  Therefore,  the  mean  scattered  pressure  measured 
experimentally  is  compared  with  the  reflection  coefficients  calculated  with  the  Somigliana 
boundary  integral  technique,  as  this  numerical  approach  includes  all  multiple  scattering  and 
wave  conversions  at  the  interface. 

Case:  A  =  0.71a 

Figure  7  shows  one  realization  of  the  interface  given  an  incident  pulse  with  a  center  frequency 
of  1.5  MHz  and  a  half-power  width  of  0.25MHz.  This  realization  corresponds  to  a  beam 
impinging  on  the  surface  with  an  incident  angle  of  20°.  The  source  pulse  along  with  the 
pulse  reflected  from  a  plane  interface  are  also  shown,  in  which  case  energy  travels  only  in 
the  specular  direction.  Referring  to  the  polar  seismogram,  it  is  clear  that  the  large  scale 
surface  roughness  scatters  energy  over  most  forward  and  back  scattering  angles.  The  energy 
is  spread  over  a  large  time  interval  and  amplitudes  vary  rapidly  as  a  function  of  scattering 
angle.  In  general,  it  is  difficult,  given  this  single  model,  to  determine  quantitatively  which 
scattering  mechanisms  are  working  at  the  surface. 

Our  first  continuous  wave  analysis  is  carried  out  at  1.5  MHz,  the  center  frequency  of  the 
seismogram  above,  so  that  A  =  .71a  =  1.00  mm.  If  one  were  to  look  at  a  single  experimental 
realization  of  the  fluid-glass  surface  at  each  of  four  incident  beam  angles:  0°,  20°,  30°,  and 
60°  one  would  find  that  the  amplitudes  in  each  realization  vary  strongly  as  a  function  of 
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scattering  angle,  6s.  Averaging  over  a  finite  number  of  these  realizations  the  mean  reflection 
coefficient  is  obtained.  The  total  mean  reflection  coefficients  for  both  a  Gaussian  and  an 
exponential  surface  are  given  in  Figures  8  through  11. At  the  bottom  of  these  figures,  the 
experimental  mean  reflection  coefficients  are  given  along  with  the  SD  of  the  finite  average, 
showing  the  deviation  of  these  reflection  coefficients  from  the  final  mean  reflection  coefficient, 
which  would  correspond  to  a  full  ensemble  of  realizations.  Negative  scattering  angles  (0,  <0) 
correspond  to  backscattering  in  all  plots.  We  also  stress  that  given  the  incident  wavelength, 
the  surface  is  extremely  irregular  and  the  specular  reflection  is  largely  disrupted. 

Figure  8  shows  the  total  mean  scattered  pressure  as  a  function  of  scattering  angle  given 
a  normally  incident  acoustic  beam.  Upon  comparing  the  numerical  and  experimental  data 
it  is  clear  that  the  2-D  numerical  results  for  a  Gaussian  interface  can  give  insight  into  the 
3-D  experimental  results.  The  fluctuations  in  the  data  are  mostly  within  one  standard 
deviation  of  the  finite  average.  There  are  a  number  of  interesting  aspects  in  the  curves.  The 
experimental  data  shows  a  strong  peak  amplitude  propagating  back  towards  the  source  at 
=  0°,  and  this  is  predicted  by  the  numerical  reflection  coefficient.  The  width  of  this  peak 
is  approximately  35°.  Notice  that  there  is  also  considerable  scattering  at  all  angles,  dropping 
oflF  linearly  with  increasing  scattering  angle.  The  exponential  surface  also  predicts  portions 
of  the  amplitude  distribution,  but  its  higher  frequency  component  appears  to  destroy  the 
enhancement  of  energy  scattered  back  towards  the  source. 

A  similar  form  of  scattering  is  exhibited  in  Figures  9  and  10,  which  show  the  mean 


19 


scattered  pressure  for  an  incident  angle  of  20°  and  30°,  respectively.  Again  the  numerical 
results  for  the  Gaussian  interface  are  similar  in  nature  to  the  experimental  results.  There  are 
also  some  remarkable  aspects  to  these  two  curves.  First,  both  reflection  coefficients  contain 
two  dominant  peaks.  One  peak  is  broad  and  occurs  in  the  forward  scattered  direction 
while  the  other  peak  is  much  narrower  and  occurs  in  the  retroreflective  direction,  Og  =  -^o- 
The  experimental  results  show  a  consistent,  yet  less  distinct,  peak  in  the  retroreflective 
direction.  This  peak  clearly  loses  amplitude  as  the  incident  angle  is  increased,  sinking  further 
into  the  surrounding  reflection  coefficient.  A  second  aspect  is  that  both  curves  show  a 
strong  asymmetry.  However,  upon  comparing  the  curves,  the  2-D  numerical  model  shows 
more  backscattering  and  less  forward  scattering  than  the  3-D  ultrasonic  data.  This  trend 
becomes  more  prominent  as  the  incident  angle  is  increased.  The  exponential  curve  follows  the 
Gaussian  reflection  coefficient  closely,  although  it  again  shows  a  less  distinct  retroreflective 
peak. 

Figure  11  gives  the  mean  scattered  pressure  for  a  beam  incident  at  60°.  In  this  case  there 
are  no  distinct  signs  of  enhanced  backscattering  in  either  the  numerical  or  the  experimental 
data.  However,  energy  is  scattered  uniformly  over  most  backscattering  angles.  This  energy 
does  not  drop  off  until  the  retroreflective  angle  is  exceeded  in  the  backscattering  region. 
Most  dramatic  is  the  continuation  of  the  trend  observed  at  the  smaller  incident  angles  above. 
Speciflcally,  the  numerical  data  clearly  shows  a  slower  falloff  in  backscattering  and  a  slower 
rise  in  forward  scattering  than  exhibited  by  the  3-D  experimental  data. 
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Case:  A  =  2.14a 


Figure  12  shows  one  realization  of  the  interface  given  an  incident  pulse  with  a  center  frequency 
of  0.5  MHz  and  a  half-power  width  of  250  kHz.  This  realization  corresponds  to  a  beam 
impinging  on  the  surface  at  a  20°  incident  angle.  The  pulse  reflected  from  a  plane  interface 
is  shown,  with  the  energy  again  arriving  only  in  the  specular  direction.  Referring  to  the 
scattered  seismogram,  it  is  clear  that,  even  at  this  lower  frequency,  energy  is  scattered  over 
most  forward  and  backscattered  angles.  Given  this  one  deterministic  case,  it  is  again  difficult 
to  determine  quantitatively  the  scattering  mechanisms  operating  at  this  frequency. 

The  second  continuous  wave  analysis  was  carried  out  at  0.5  MHz,  the  center  frequency 
of  the  seismogram  above,  so  that  A  =  2.14a  =  3.0  mm.  Once  again,  if  one  were  to  look  at  a 
single  realization,  the  amplitudes  for  each  realization  vary  strongly  as  a  function  of  scattering 
angle,  although  not  as  strongly  as  the  A  =  .71a  case.  Figure  13  shows  the  comparison  be¬ 
tween  the  numerical  and  experimental  mean  reflection  coefficients  given  a  normally  incident 
beam.  The  2-D  numerical  results  for  a  Gaussian  interface  are  able  to  predict  portions  of 
the  experimental  observations.  In  this  case,  much  of  the  experimental  data  sits  within  one 
standard  deviation  of  the  finite-average.  Comparing  these  curves  to  the  curves  for  A  =  0.71a, 
a  number  of  distinct  differences  are  apparent.  Most  noticeable  is  the  widening  of  the  retrore- 
flective  peak  width  from  about  35°  to  greater  than  60°.  This  widening  is  apparent  in  both  the 
experimental  and  the  numerical  data.  The  reflection  coefficient  for  an  exponential  interface 
again  shows  much  lower  retroreflectance  than  for  the  Gaussian  interface. 
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Figures  14  and  15  both  show  that  the  numerical  results  over  a  Gaussian  interface  predict 
the  asymmetric  trends  in  the  experimental  data  for  incident  angles  of  20°  and  30°,  respec¬ 
tively.  However,  distinct  differences  do  occur  between  the  two  curves.  First,  as  the  incident 
angle  increases,  the  2-D  numerical  results  again  show  more  backscattering  and  less  forward 
scattering  than  the  3-D  ultrasonic  data.  A  broad  retroreflective  peak  is  both  predicted  and 
observed  at  20°  and  30°  incidence,  supporting  the  existence  of  retroreflectance.  Unfortu¬ 
nately,  the  height  of  these  peaks  are  of  the  same  order  as  the  standard  deviation  of  the 
experimental  average,  not  allowing  for  a  direct  verification  of  retroreflectance.  Numerically, 
the  exponential  interface  does  give  rise  to  a  retroreflective  peak,  although  this  peak  is  smaller 
than  the  peak  predicted  by  the  Gaussian  surface. 

Figure  16  shows  the  mean  reflection  coefficient  for  an  incident  angle  of  60°.  In  this 
case  enhanced  backscattering  is  not  predicted  numerically  or  observed  experimentally.  The 
2-D  numerical  model  again  predicts  a  slower  falloff  in  backscattering  and  a  slower  rise  in 
forward  scattering  than  exhibited  by  the  3-D  ultrasonic  data.  In  addition,  the  numerical 
model  predicts  a  much  smaller  specular  peak  than  is  observed  experimentally.  Although  the 
amplitudes  are  different,  the  numerical  curves  do  predict  the  uniform  scattering  of  energy 
into  the  fluid  above  as  shown  in  the  data. 

Retroreflectance  is  clearly  supported  by  the  ultrasonic  data  above.  However,  the  retrore¬ 
flective  peak  height  is  still  on  the  same  order  as  the  standard  deviation  of  the  finite-average 
in  each  case.  This  makes  it  difficult  to  verify  the  existence  of  “enhanced  backscattering” 
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absolutely.  For  this  reason,  data  was  recorded  near  the  retroreflective  direction  using  the 
partially  transparent  piezofilm  receiver.  The  final  average  RC  observed  with  the  piezofilm  re¬ 
ceiver  is  superposed  on  Figure  14.  The  data  has  been  scaled  to  the  amplitudes  received  with 
the  flat-bottomed  transducers.  The  scattered  pressure  was  measured  between  the  backscat- 
tering  angles  of  40°  and  5°  (—40°  <  Os  <  —5°),  and  65  surface  realizations  were  averaged.  In 
this  case,  a  distinct  peak  is  observed  in  the  retroreflective  direction  with  a  slightly  narrower 
form  than  the  numerically  generated  peak.  The  size  of  the  retroreflective  peak  is  larger  than 
the  corresponding  SD  of  the  average,  strongly  supporting  the  the  presence  of  an  enhancement 
of  retroreflected  energy  due  to  multiple  scattering  from  the  glass  interface. 

Upon  a  careful  review  of  the  numerical  and  experimental  curves  above  there  are  finer 
features  in  the  experimental  data  which  may  be  of  importance.  For  instance,  the  experi¬ 
mental  observations  show  a  distinct  change  in  the  variance  of  the  amplitude  measurements. 
In  the  both  the  A  =  2.14a  and  the  A  =  .71a  case,  the  variance  of  the  mean  reflection  coeffi¬ 
cient  decreases  dramatically  outside  the  —60°  <  Og  <  +60°  range.  This  change  in  variance 
does  not  appear  to  depend  strongly  on  the  incident  angle  or  the  incident  wavelength,  at 
least  in  the  frequency  range  studied.  The  change  occurs  in  conjunction  with  a  sharp  break 
and  a  distinct  drop  in  the  reflected  amplitudes  at  wider  scattering  angles,  suggesting  that 
it  may  be  associated  with  shadowing  effects  at  the  larger  scattering  angles.  Numerically, 
by  averaging  a  similar  number  of  2-D  numerical  realizations,  this  break  in  the  variance  for 
a  small  finite  average  can  be  predicted.  As  more  realizations  are  averaged,  the  sharpness 


23 


of  the  break  in  reflected  amplitudes  becomes  smoothed.  Surprisingly,  the  model  with  an 
exponential  correlation  function  predicts  the  point  of  this  break  better  than  the  Gaussian 
model. 

Upon  comparing  the  experimental  and  numerical  curves,  additional  discrepancies  can  be 
identifled  between  the  numerically  predicted  and  the  experimentally  observed  amplitudes. 
The  most  prominent  discrepancy  is  the  appearance  of  a  large  amplitude  specular  peak  at 
small  incident  angles  in  the  3-D  data.  While  a  peak  is  clearly  present  in  the  specular  direction 
by  a  20°  incident  angle  in  the  each  3-D  case,  a  distinct  specular  peak  is  only  present  for 
9q  =  30°  in  the  A  =  2.14a  numerical  case.  Discrepancies  such  as  this  can  be  useful  for 
understanding  the  effect  of  2-D  and  3-D  structures  on  seismic  motions. 

Discussion  and  Conclusions 

In  this  study,  we  were  able  to  generate,  within  reasonable  accuracy,  a  3-D  characterized  ran¬ 
dom  interface.  An  interface  with  approximately  a  Gaussian  surface  height  distribution  and  a 
Gaussian  correlation  function  was  generated  using  a  glass  etching  procedure  and  photoresist 
templates.  The  resulting  surface  distribution  was  confirmed  using  surface  profilometer  mea¬ 
surements.  Scattered  pressures  were  then  acquired  over  this  surface  and  compared  directly 
to  numerical  results  for  a  2-D  interface  with  the  same  statistical  parameters. 

Generally,  the  numerically  derived  mean  reflection  coefficients  calculated  over  an  acoustic- 
elastic  interface  show  retroreflective  trends  similar  to  those  observed  for  the  SH  and  P-SV 
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cases.  First,  the  width  of  the  retroreflective  peak  appears  to  be  consistent  with  the  multiple 
scattering,  constructive  phase  argument  summarized  in  the  introduction.  In  this  case,  when 
the  wavelength  is  increased,  both  the  numerically  derived  curve  and  the  experimental  data 
show  an  increase  in  peak  width,  from  35°  at  A  =  .71a,  to  greater  than  60°  at  A  =  2.14a. 
Second,  as  the  incident  angle  is  increased,  the  retroreflective  peak  amplitude  tends  to  decrease 
relative  to  the  remaining  portion  of  the  reflection  coefficient.  The  retroreflective  peak  in 
both  the  experimentally  and  numerically  derived  curves  disappears  at  an  incident  angle 
approximately  equal  to  the  30°  rms  slope  of  the  interface.  Finally,  although  not  studied 
directly,  it  seems  likely,  based  on  the  work  of  Schultz  (1993),  that  the  retroreflective  peak 
height  will  tend  to  decrease  as  the  impedance  contrast  is  lowered  and  more  energy  is  allowed 
to  penetrate  the  interface.  Along  the  same  lines,  the  retroreflective  peak  amplitude  is  likely 
to  diminish  as  the  rms  slope  of  the  interface  is  decreased,  since  not  as  many  time-reversed 
paths  can  be  obtained  with  the  lower  slopes. 

Near  normal  incidence,  we  have  shown  that  the  2-D  numerical  results  for  scattering  from  a 
randomly  irregular  interface,  with  Gaussian  statistics,  can  give  insight  into  the  3-D  scattering 
observed  experimentally.  At  larger  incident  angles,  the  2-D  numerical  curves  deviated  from 
the  experimental  curves  significantly.  More  specifically,  as  the  incident  angle  was  increased  to 
greater  than  20°,  the  2-D  models  predicted  a  slower  decrease  in  backscattering  and  a  slower 
increase  in  forward  scattering  than  observed  in  the  3-D  experimental  data.  The  trend  became 
dramatic  as  the  incident  angle  was  increased  to  60°.  Since  the  surface  in  this  experiment  was 
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well  characterized,  this  appears  to  be  the  result  of  inherent  differences  between  2-D  and  3-D 
scattering  mechanisms.  More  precisely,  as  the  incident  angle  is  increased  the  extra  degree 
of  spatial  freedom  a.ssociated  with  a  3-D  interface  appears  to  favor  the  propagation  of  more 
forward  scattered  energy  and  less  backscattered  energy  relative  to  a  2-D  interface.  Although 
the  falloff  in  amplitudes  was  different  between  the  2-D  and  3-D  interfaces,  the  2-D  numerical 
results  did  predict  the  general  trend  of  the  experimental  scattering.  Both  results  showed  that 
at  large  incident  angles  energy  is  scattered  uniformly  over  most  scattering  angles,  resulting 
in  negative  phase  velocities,  large  phase  velocities,  and  a  large  amount  of  interference  in 
seismic  data  recorded  in  the  fluid  above. 

Although  tests  were  not  carried  out  on  epoxy  surfaces  in  this  study,  an  epoxy  surface 
was  generated  for  the  purpose  of  profilometer  measurements,  by  using  the  irregular  glass 
surface  as  a  mold.  Future  work  may,  therefore,  include  studies  to  determine  how  material 
properties  affect  scattering  from  surfaces  with  identical  height  distributions.  Also,  a  similar 
etching  process  may  be  used  to  study  interfaces  with  differing  statistics. 
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Figure  1:  Peak-valley  sequence  along  an  interface  showing  an  example  of  the  time- reversed 
paths  which  lead  to  enhanced  backscattering.  The  solid  line  shows  a  forward  scattered 
path  while  the  dashed  line  shows  the  corresponding  time-reversed  path.  These  two  paths 
interfere  constructively  to  give  an  increase  in  amplitude  diflPracted  back  in  the  direction 
of  the  source. 


Figure  2.  The  six  templates  used  to  generate  the  random  surface  used  in  these  experiments. 
Each  circle  has  a  19.5  cm  diameter  to  match  the  glass  surface.  Each  template  corresponds 
to  one  standard  deviation  of  depth  and  each  template  was  exposed  to  high  velocity 
particles  for  the  same  amount  of  time.  The  numbering  shows  the  order  in  which  the 
templates  were  applied. 


SIDE  VIEW 


Figure  3:  The  geometry  used  to  experimentally  measure  the  scattering  properties  of  a  given 
random  surface.  The  source  is  held  stationary  at  one  incident  angle  while  the  detector 
is  stepped  in  a  semi-circle  about  the  random  surface.  This  then  gives  one  realization  of 
that  surface. 
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Figure  4:  Histogram  plotting  surface  height  data.  The  target  surface  height  distribution 
(a)  is  Gaussian  with  a  standard  deviation  1  mm.  The  surface  height  distribution  (b) 
based  on  profilometer  measurements  (squares)  of  the  glass  surface  is  shown  along  with  a 
best  fitting  Gaussian  curve  (solid  line)  which  has  a  standard  deviation  of  0.6  mm.  This 
histogram  was  plotted  using  40000  surface  profilometer  measurements. 
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Figure  5:  The  interface  autocorrelation  function.  The  target  autocorrelation  function  (a) 
is  a  Gaussian  function  with  a  correlation  length  e~^  of  1.0  mm.  The  actual  autocor¬ 
relation  function  of  the  glass  block  (b)  as  calculated  from  profilometer  measurements 
has  a  correlation  length  of  1.4  mm.  The  surface  profilometer  measurements  (solid  line) 
are  compared  with  Gaussian  (crosses)  and  exponential  (circles)  functions  with  the  same 
correlation  lengths  of  1.4  mm. 
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Figure  6.  The  surface  height  distribution  based  on  profilometer  measurements,  (a)  gives  a 
grayscale  plot  of  the  surface,  (b)  gives  a  profile  across  the  surface  as  marked  in  (a),  and 
(c)  shows  a  numerically  generated  surface  with  the  statistics  given  in  Figures  4  and  5. 
Both  Gaussian  (solid  line)  and  exponential  (dashed  line)  correlation  functions  are  shown. 
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Figure  7:  Experimental  seismogram  recorded  over  the  glass  model,  (a)  shows  the  seismic 
data  recorded  as  a  function  of  angle  over  the  irregular  glass  surface  with  A  =  .71a  given 
an  acoustic  beam  incident  at  20°.  The  arrow  shows  the  retroreflective  direction.  The 
source  wavelet  (b)  and  the  specular  reflection  (c)  recorded  over  a  plane  interface  are  also 
shown.  In  the  plane  layer  case  the  only  observable  energy  is  in  the  specular  direction. 


Figure  8:  The  2-D  mean  reflection  coefficient  (a)  calculated  numerically  over  the  Gaussian 
(solid  line)  and  exponential  (dashed  line)  surfaces  given  a  normally  incident  source  beam 
with  A  =  .71a.  The  3-D  reflection  coefficient  (b)  recorded  over  the  experimental  interface 
(circles)  and  the  standard  deviation  of  the  finite  average  (plus)  are  also  shown.  The 
experimental  results  correspond  to  30  surface  realizations. 


Figure  9:  Similar  to  Figure  8,  except  the  incident  angle  is  now  20°  and  results  correspond 
to  20  realizations. 
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Figure  10;  Similar  to  Figure  8,  except  the  incident  angle  is  now  30°  and  results  correspond 
to  10  realizations. 
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Figure  11:  Similar  to  Figure  8,  except  the  incident  angle  is  now  60°  and  results  correspond 
to  10  realizations. 
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Figure  12:  Similar  to  Figure  7,  except  A  =  2.14a. 
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Figure  13.  The  2-D  reflection  coefficient  (a)  calculated  numerically  over  the  Gaussian  (solid 
line)  and  exponential  (dashed  line)  surfaces  given  a  normally  incident  source  beam  with 
^  2.14a.  6s  is  the  angle  of  forward  (positive)  and  back  (negative)  scattering.  The 

3-D  reflection  coefficient  (b)  recorded  over  the  experimental  interface  (circles)  and  the 
standard  deviation  of  the  flnite  average  (plus)  are  also  shown.  The  experimental  results 
correspond  to  30  surface  realizations. 
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Figure  14:  Similar  to  Figure  13,  except  the  incident  angle  is  now  20°  and  results  correspond 
to  20  realizations. 


Figure  15:  Similar  to  Figure  13,  except  the  incident  angle  is  now  30°  and  results  correspond 
to  10  realizations. 
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Figure  16:  Similar  to  Figure  13,  except  the  incident  angle  is  now  60°  and  results  correspond 
to  10  realizations. 
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CRUSTAL  REFLECTIONS  AND  THE  NATURE  OF 

REGIONAL  P  CODA 


Summary 

F-K  spectra  for  regional  P  coda  are  presented  for  events  recorded  at  the  Scandinavian 
NORESS,  FINESA,  and  ARCESS  arrays  and  the  New  England  NYNEX  array.  The  spectra 
are  dominated  by  on-azimuth  energy  with  apparent  velocities  between  Pn  (or  faster)  and  Lg. 
Following  this  analysis,  reflection  coefficients  calculated  with  an  efficient  boundary  integral 
scheme  are  used  to  study  the  role  irregular  interfaces  play  in  the  creation  of  regional  P  coda. 
We  find  that  observed  crustal  scattering  in  these  regions  is  strikingly  consistent  with  P-P 
and  P-SV  scattering  from  the  2-D  irregular  Moho  and  even  more  consistent  with  scattering 
from  a  2-D  irregular  near  surface  interface. 

Introduction 

In  recent  years,  understanding  regional  seismograms  has  attracted  considerable  interest,  in 
part  because  of  a  desire  to  monitor  small  underground  nuclear  tests  at  regional  distances. 
One  of  the  characteristics  of  regional  seismograms  is  coda,  a  train  of  waves  without  evident 
phases  following  the  P  and  S  arrivals.  Previous  work  has  shown  that  both  regional  P  coda 
and  S  coda  recorded  in  the  crust  are  more  complicated  than  the  coda  predicted  for  a  plane 
layered  crust.  Dainty  and  Toksdz  (1990)  applied  F-K  analysis  at  Scandinavian  arrays.  They 
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found  that  at  later  times  in  the  S  coda,  energy  tends  to  arrive  from  all  directions.  This 
energy  is  consistent  with  the  isotropic  “standard  model”  of  Aki  and  Chouet  (1975),  modified 
for  the  regional  case  by  considering  the  scattering  to  be  Lg  to  Lg  (e.g.,  Herrmann,  1980). 
Teleseismic  P  coda  appears  to  be  similar  in  many  respects  to  S  coda  (Dainty,  1990;  Dainty 
and  Harris,  1989;  Dainty  and  Toksoz,  1990).  Regional  P  coda,  on  the  other  hand,  is  more 
enigmatic  and  one  model  does  not  appear  to  explain  the  observations.  The  final  model  of 
regional  P  coda  is  likely  a  combination  of  various  crustal  scattering  mechanisms. 

Scattering  due  to  irregular  interfaces  may  give  a  better  understanding  of  the  incoherent, 
low  amplitude  arrivals  which  tend  to  compficate  these  seismograms.  Recently,  boundary 
integral  techniques  have  been  shown  to  be  an  effective  tool  in  identifying  which  scattering 
mechanisms  exist  at  a  given  interface.  These  efficient  formulations  have  allowed  the  devel¬ 
opment  of  fast  algorithms  which  can  model  statistical  reflection  coefficients  over  randomly 
irregular  interfaces  (Schultz  and  Toksoz,  1993;  Schultz  and  Toksoz,  1994). 

In  this  study,  we  illustrate  how  the  reflection  coefficients  developed  by  Schultz  and  Toksoz 
(1994)  can  be  used  to  gain  insight  into  crustal  scattering  mechanisms.  We  first  describe  the 
general  nature  of  regional  P  coda  as  recorded  at  Scandinavian  (Mykkeltveit  et  a/.,  1990) 
and  New  England  (Dainty  and  Battis,  1989)  arrays  in  terms  of  azimuth  and  phase  velocity 
(angle  of  incidence)  of  energy.  We  characterize  P  coda  using  an  F-K  spectral  approach. 
We  then  examine  numerically  derived  reflection  coefficients  from  Schultz  and  Toksoz  (1994) 
to  determine  whether  scattering  from  interface  irregularities  is  consistent  with  the  observed 
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F-K  spectra. 


The  Nature  of  Regional  P  Coda 

Regional  P  coda  consists  of  a  train  of  waves  after  Pn  and  Pg  and  before  Sn  and  Lg  at 
regional  distances  (100-2000  km).  On  the  instruments  used  in  this  study,  frequencies  are 
in  the  range  1-20  Hz.  Explanations  for  P  coda  include  the  scattering  of  P  and  S  waves 
from  distributed  small  heterogeneities  in  a  similar  manner  to  S  coda,  and  multiple  crustal 
reflections  predicted  by  synthetic  seismogram  calculations  in  laterally  homogeneous  layered 
media.  The  results  of  Dainty  and  Toksoz  (1990)  show  that  P  coda  cannot  be  explained  by 
the  same  model  as  S  coda  and  a  composite  model  which  includes  various  forms  of  scattering 
may  be  required  to  explain  P  coda  observed  at  arrays  in  Scandinavia.  Such  a  composite 
model  must  consider  the  possibility  of  scattering  from  interface  irregularities.  Therefore,  in 
this  section,  array  data  is  studied  in  more  detail  to  determine  whether  scattering  from  an 
irregular  Moho  and  near  surface  layers  is  consistent  with  the  observed  P  coda.  We  turn  our 
attention  to  the  analysis  of  regional  P  coda  (100-2000  km)  recorded  at  the  Scandinavian 
NORESS,  FINESA,  and  ARCESS  arrays  and  the  New  England  NYNEX  array. 

Figure  1  gives  a  synthetic  vertical  seismogram  section  for  a  plane  layered  crust  recorded 
along  an  array  of  receivers  at  distances  of  100  to  500  km  from  a  vertical  point  force  applied  at 
a  depth  of  .04  km  (from  Toksoz  et  al,  1990).  The  model  summarized  in  Table  1  is  based  on 
the  NORESS  crustal  model  of  Sellevol  and  Warwick  (1971)  where  Q  =  Qp^sf^’,  and  serves  to 
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Table  1:  NORESS  Crustal  Model 


Depth  (km) 

Vj,  (km/s) 

Vs  (km/s) 

P  (gm/cc) 

Qp 

Qs 

C 

0.0 

5.2 

3.2 

2.6 

50 

25 

0.5 

3.0 

6.0 

3.5 

2.7 

200 

100 

0.8 

17.0 

6.51 

3.8 

2.9 

1000 

500 

0.0 

38.0 

8.05 

4.67 

3.3 

2000 

1000 

0.0 

illustrate  the  effects  discussed.  Referring  to  Figure  1  the  synthetic  results  show  that  a  plane 
layered  crust  results  in  P  coda  energy  which  will  arrive  on  azimuth,  with  the  predominant 
energy  showing  phase  velocities  between  Pg  and  Pn.  P-SV  bounces  are  included  but  they 
have  high,  Pg-Pn,  phase  velocities.  This  behavior  is  caused  by  rays  with  multiple  bounce 
points  in  the  crust,  mainly  on  the  Moho  and  the  free  surface  where  at  least  one  leg  is  P.  At 
larger  offsets  the  energy  in  P  coda  converges  to  a  Pg  phase  velocity.  Energy  is  not  expected 
at  velocities  lower  than  Pg  or  higher  than  Pn.  Also,  the  P  coda  energy  in  this  plane  layer 
model  is  coherent  with  arrival  groups  that  may  be  correlated  across  traces.  This  coherence  is 
not  observed  in  real  seismograms.  The  phenomenon  of  multiple  post-critical  crustal  bounces 
potentially  distinguishes  regional  P  coda  from  teleseismic  P  coda,  where  crustal  reflections 
will  be  precritical  and  of  low  energy. 

F-K  Analysis 

F-K  analysis  is  used  to  find  the  velocity  and  azimuth  of  the  dominant  and  secondary  energy 
observed  in  the  P  coda  at  the  four  arrays;  NORESS,  ARCESS,  (Mykkeltveit  et  aL,  1990) 
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Table  2:  Event  parameters  and  measured  azimuths  and  velocities  from  frequency- 
wavenumber  analysis 


mBSSmi 

IBH 

P  Az 

1st  Peak 

2nd  Peak 

3rd  Peak  I 

Az 

BanasxoB 

Si 

BitasrxiBi 

Az 

KfiWaig 

HBHBI 

7 

420 

3.8 

8.0 

BilBi! 

■IB^^B 

M 

370 

2.4 

■TCITM 

■Xik&liB 

8.5 

E 

1160 

2.2 

^ZSiB 

BTVff'HB 

4.5 

N 

2.0 

■IQSiB 

BElS]Bi 

5.8 

■j^ZuiB 

HKE^^H 

BBHBH 

M 

420 

4.0 

7.9 

■QEBi 

HHHHI 

M 

430 

4.0 

6.2 

BOOM 

HHHH 

M 

240 

4.0 

HBBHB 

M 

230 

4.0 

BiliZilB 

7.25 

HBIB^I 

hesembb 

M 

220 

4.0 

■EliSiB 

IBBBB 

11:06:26 

M 

340 

4.0 

K£SiB 

IHHHi 

HIHBB 

BLUifisisH 

M 

300 

4.0 

6.0 

KESB 

HHBHi 

HBiBHI 

■HSIiloS) 

M 

320 

2.4 

6.0 

E 

360 

6.5 

kiebi 

7.6 

Mill  f 

■eeeesxih 

E 

5.0 

II^SQiB 

7.5 

KE&ZaiB 

M 

1070 

3.0 

8.5 

^HIBI 

BBBBBl 

hzbemh 

M 

390 

3.0 

7.8 

BEBiB 

6.25 

HHBEBB 

NYNEX 

R 

165 

15.0 

91. O*' 

■ZZliliB 

15.7 

40.0° 

R 

198 

HHflQHI 

tEMSiM 

KUiZiliB 

HHZBHI 

li^^BBI 

HHHHH 

R 

189 

12.0 

■QQQiB 

IHBBI 

'BHIB 

IHBIH^I 

mamm 

159 

7.0 

BTflTliB 

■ipyiM 

IHHQQHI 

IBBBB 

■^■HB 

240 

8.0 

;E£SQiB 

IBKBBH 

IBBBBi 

IBBBi 

IBBBBB' 

“Calculated  from  S-P  time,  or  taken  from  Helsinki  Bulletin,  or  taken  from  shot  table. 

*’E=  Earthquake,  M  =  Mine  Blast,  N  =  Presumed  Nuclear  Test,  R  =  Refraction  Shot, 
?  =  Unknown 

FINES  A,  (Uski,  1990)  and  NYNEX  (Battis,  1990).  A  combined  total  of  21  events  were 
studied  consisting  of  a  combination  of  mine  blasts,  nuclear  blasts,  refraction  shots  (NYNEX) 
and  earthquakes  (Table  2).  Figure  1  gives  an  example  of  the  processing  and  analysis,  where 
the  seismogram  corresponds  to  the  vertical  displacement  recorded  at  the  center  station  of 
NORESS%.  The  event  was  a  quarry  blast  at  Titania  Quarry,  which  is  393  km  away  from 
the  array  with  an  azimuth  of  228°.  The  algorithm  used  for  F-K  analysis  is  that  described 
in  Dainty  and  Toksoz  (1990).  The  seismogram  is  filtered  in  a  narrow  pass-band,  -f/  —  10%, 
around  a  central  frequency  of  3  Hz  and  the  frequency-wavenumber  spectrum  is  calculated 
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using  the  Maximum  Likelihood  Method  (MLM)  (Capon,  1969).  A  window  of  20  sec  duration 
was  used  for  NORESS,  ARCESS,  and  FINES  A  data,  10  sec  for  NYNEX  data;  the  windows 
were  picked  by  hand.  While  the  analysis  method  allows  overlapping  windows,  this  was  not 
used  in  this  study.  Frequencies  were  chosen  to  be  within  the  range  that  had  significant 
power  in  the  data,  usually  the  peak  frequency.  The  MLM  method  gives  the  resolution  in 
wavenumber  needed  for  this  study;  a  narrow  band  was  used  because  our  experience  with 
broadband  MLM  is  that  it  is  very  computer  intensive  and  narrow  band  MLM  was  used 
with  success  by  Dainty  and  Toksdz  (1990).  AU  velocities  and  azimuths  reported  here  were 
determined  by  this  method. 

The  results  from  the  F-K  analysis  showed  peaks  at  close  to  the  source  azimuth.  This 
encouraged  us  to  use  semblance  analysis  to  determine  the  apparent  velocity  of  P  coda  energy 
as  a  function  of  time  along  the  source  azimuth  (228°  in  Figure  2)  (Dainty  and  Schultz,  1993). 
The  semblance  S  at  slowness  s  was  calculated  by  the  standard  formula  (Kanasewich,  1975, 
p.  252) 

<7™ 

where  fij  are  the  seismogram  values  at  array  station  j  and  time  index  i,  bi  are  the  values 
of  the  beam  at  slowness  s  (and  the  source  azimuth)  at  time  index  i,  M  is  the  total  number 
of  array  stations,  and  the  summation  over  i  represents  the  time  window  over  which  the 
semblance  is  calculated.  This  quantity  was  calculated  in  2  sec  windows  overlapped  by  1  sec 
and  a  slowness  range  of  0.0  to  0.1  sec/km  at  the  source  azimuth.  The  resulting  grid  of  values 
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was  contoured  to  show  where  coherent  energy  was  arriving  in  time  (cf.  Figure  2).  Semblance 
has  a  range  of  zero  to  one  and  is  a  very  sensitive  indicator  of  the  presence  of  coherent  energy; 
in  Figure  2  light  areas  in  the  semblance  plot  have  a  semblance  of  0.4  and  greater.  However,  it 
must  be  noted  that  velocities  correponding  to  peaks  on  the  semblance  plot  shown  in  Figure  2 
will  only  be  correct  if  the  energy  is  arriving  at  exactly  the  azimuth  used  to  construct  the 
beams  (228°).  For  this  reason  velocities  were  not  picked  from  the  semblance  plots. 

F-K  analysis  of  P  coda  energy  assists  in  the  elucidation  of  crustal  scattering  mechanisms. 
This  preliminary  analysis  involved  picking  by  hand  the  phase  velocities  and  azimuthal  devi¬ 
ations  of  energy  peaks  in  the  F-K  spectra.  The  peaks  were  chosen  based  on  several  criteria. 
Generally,  one  or  two  peaks  were  picked  on  a  plot  and  for  one  event  (Figure  2)  three  peaks 
were  picked.  The  peaks  selected  were  the  highest  peaks,  at  least  one  contour  level  above  the 
surrounding  noise  with  phase  velocities  faster  than  2  km/s.  The  largest  peak  was  chosen  as 
a  primary  peak  and  the  remaining  peak(s)  were  labeled  as  secondary.  No  secondary  peak 
was  more  than  6  db  below  the  primary,  and  all  but  one  were  within  3db.  Figure  3  gives  the 
histograms  for  the  phase  velocity  and  azimuthal  distribution  of  the  data  and  includes  both 
primary  and  secondary  peaks.  The  peaks  cluster  around  the  source  azimuth  (most  peaks 
within  20°).  This  contrasts  with  the  S  coda  where  the  energy  is  spread  over  all  azimuths. 
This  F-K  analysis  shows  energy  arriving  between  Pg  and  Pn  phase  velocities  as  predicted 
by  the  plane  layered  geometry.  Peak  energy  is  observed  at  Pn  phase  velocity  and  greater. 
Not  predicted  by  the  plane  layer  model,  a  large  amount  of  energy  arrives  with  slower  phase 
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velocities,  between  Lg  and  Pg,  with  primary  peak  energy  arriving  at  Sn.  A  peak  number 
of  secondary  arrivals  are  also  observed  with  approximate  Lg  velocity.  Each  of  these  slower 
arrivals  suggests  that  P  waves  are  converted  during  some  leg  of  the  path  into  post  critical- 
shear  energy  which  then  becomes  favored  for  lateral  propagation.  The  semblance  analysis 
generally  confirms  the  F-K  observations  as  Pn  and  Pg  energy  tends  to  arrive  early  in  the 
P  coda  while  Lg  energy  tends  to  arrive  late.  This  time  delay  is  expected  with  conversion 
to  shear  along  at  least  one  leg  of  the  path.  In  general  this  F-K  analysis  shows  that  P  coda 
energy  is  usually  dominated  by  a  few  peaks  which  arrive  with  different  phase  velocities  of  ap¬ 
proximately  the  same  azimuth  rather  than  a  single  phase  velocity  with  substantially  different 
azimuths,  characteristic  of  S  coda. 

Although  volume  scattering  appears  to  explain  the  scattering  responsible  for  S  coda,  the 
wide  range  of  phases  and  the  absence  of  diffuse  scattering  over  a  wide  range  of  azimuths  in 
P  coda  suggests  a  different  scattering  mechanism.  Many  mechanisms  may  be  responsible  for 
the  observation  of  slower  phases  in  P  coda  at  these  arrays.  Sanchez-Sesma  and  Campillo 
(1991),  among  others,  have  shown  theoretically  that  even  simple  topographic  features  near  an 
array  can  cause  great  variability  in  recorded  amplitudes  as  a  function  of  frequency,  incident 
angle,  and  receiver  location  with  respect  to  the  topography  and  observations  confirm  this 
(Bannister  et  al,  1990,  Gupta  et  di,  1990).  Heterogeneities  near  the  receiver  array  can 
also  play  havoc  with  any  F-K  interpretation  as  near  field  interference  spreads  and  distorts 
coherent  arrivals  on  seismograms.  Another  mechanism  which  also  must  be  considered  is 
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wave  conversion  at  irregular  crustal  interfaces  such  as  the  Moho  and  shallow  near  surface 
interfaces. 

We  assume  that  F-K  analysis  emphasizes  energy  consisting  of  plane  wave  components, 
and  thus  accentuates  waves  scattered  more  than  a  few  wavelengths  from  the  array.  We  now 
show  that  the  energy  scattered  at  irregular  crustal  boundaries  is  consistent  with  the  energy 
observed  at  the  NORESS,  FINESA,  ARCESS,  and  NYNEX  arrays. 

The  Role  of  an  Irregular  Moho 

We  will  analyze  the  F-K  results  in  terms  of  scattering  of  plane  waves  incident  upon  an 
irregular  interface,  taken  here  to  be  the  Moho.  The  comparison  between  reflection  coefficients 
and  F-K  spectra  is  based  on  the  relation  Vapp  =  Vtrue/sinO,  where  Vapp  is  the  apparent  (phase) 
velocity  measured  across  the  array,  Vtme  is  the  velocity  of  the  surface  or  crustal  material, 
and  9  is  the  angle  of  propagation  to  the  vertical.  We  utilize  the  Somigliana  approach  of 
Schultz  and  Toksoz  (1993,  1994)  to  calculate  the  effects  of  scattering.  This  approach  models 
all  waves  scattered  at  an  interface,  including  all  converted,  multiple  scattered,  and  interface 
waves  generated  along  the  interface.  The  scattered  energy  at  any  point  along  the  interface  is 
presented  in  the  form  of  amplitude  reflection  coefficients,  showing  on  average  the  amplitude 
distribution  for  waves  propagating  away  from  the  irregular  interface. 

Characterizing  the  interface  with  specific  Gaussian  statistics  the  interface  roughness  can 
be  expressed  in  terms  of  the  standard  deviation  of  interface  height,  8,  and  an  autocorrelation 
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length,  a,  which  represents  approximately  the  average  distance  between  adjacent  peaks  and 
valleys  along  the  interface.  The  average  rms  slope  of  the  interface  can  be  expressed  as 

(j)  =  (1) 

a 

We  consider  the  numerical  solutions  for  two  specific  crustal  boundaries  which  may  be  highly 
irregular.  The  first  case  is  an  irregular  Moho  where  we  assume  that  the  lower  crust  has  a 
P  wave  velocity  of  6700  m/s,  an  SV  velocity  of  3900  m/s,  and  a  density  of  3.0.  The  mantle 
is  given  a  P  wave  velocity  of  8200  m/s,  an  SV  wave  velocity  of  4700  m/s,  and  a  density 
of  3.3.  In  these  cases  we  have  selected  rather  extreme  velocity  contrasts  so  that  various 
scattering  phenomena  may  be  accentuated.  The  second  case  is  a  shallow  near  surface  soil- 
basement  interface  where  the  overlying  sediment  is  given  a  P  wave  velocity  of  2000  m/s,  an 
SV  wave  velocity  of  1200  m/s,  and  a  density  of  2.0  which  represents  an  unconsolidated  or 
semiconsolidated  sediment.  The  basement  rock  is  given  a  P  wave  velocity  of  6400  m/s,  an 
SV  wave  velocity  of  3200  m/s,  and  a  density  of  2.7.  Figures  5  and  6  show  the  statistical 
reflection  coefficients  for  a  P  wave  incident  on  a  Moho  with  both  a  10°  rms  slope  and  a 
more  irregular  30°  rms  slope,  respectively.  These  reflection  coefficients  correspond  to  the 
CELSG  CL  ^  A  presented  in  Schultz  and  Toksoz  (1994),  where  both  reflected  P  and  converted  S 
amphtudes  are  plotted  at  specific  incident  angles,  ranging  from  normal  incidence  (^o  =  0°) 
to  near  grazing  angles  (^o  =  60°).  In  each  case  the  total  displacement  amplitude  is  plotted, 
including  specular  contributions.  Note  also  that  in  each  figure,  negative  scattering  angles, 
9s,  corresponds  to  backscattering. 
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A  simplified  crustal  model,  where  a  constant  velocity  crust  is  bounded  by  a  plane  free 
surface  and  an  irregular  Moho  as  shown  in  Figure  4a,  can  help  elucidate  the  generation  of  P 
coda.  In  the  case  shown,  an  incident  P  wave  encounters  an  irregular  Moho,  generating  both 
P  and  converted  S  energy  over  a  broad  range  of  directions,  including  the  specular  direction. 
Referring  to  the  reflection  coefficients  for  an  irregular  Moho,  shown  in  Figures  5  and  6,  it 
is  clear  that  in  the  3-4  Hz  frequency  range  characteristic  of  the  Scandinavian  events,  Moho 
roughness  of  100-500  m  over  1-2  km  range  can  explain  the  observed  scattering  trends.  At  the 
higher  frequencies,  typical  of  the  NYNEX  data,  30-200  m  variations  over  300-600  m  would 
be  appropriate.  P-to-P  scattering  at  most  incident  angles  favors  a  lobe  of  scattering  with  the 
peak  amount  of  energy  scattered  at  Pn  phase  velocity  (critical  angle)  and  a  full  spectrum 
of  energy  scattered  between  Pn  and  Pg  phase  velocities  (post-critical  angles)  given  most 
incident  angles.  These  waves  add  tO  energy  arriving  in  the  specular  direction,  contributing 
to  the  wavefield  with  Pn  and  Pg  phase  velocity.  At  all  incident  angles,  P-to-S  reflection 
coefficients  show  a  peak  conversion  to  phase  velocities  faster  than  Sn.  At  incident  angles 
greater  than  40°,  energy  is  converted  to  post-critical  shear  waves.  Since  the  peak  S  energy, 
Va>  Sn,  leaks  out  of  the  crust,  the  most  efficient  modes  of  propagation  are  waves  travelling 
with  Sn  phase  velocity.  Lg  velocities  are  also  favored,  although  less  energy  is  converted  to 
these  waves.  This  distribution  of  energy  corresponds  well  to  the  observed  data  where  Sn 
and  Lg  velocities  are  seen.  Moving  to  the  more  irregular  case  of  a  30°  rms  slope,  shown  in 
Figure  6,  this  post-critical  scattering  becomes  more  pronounced,  demonstrating  that  larger 
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interface  roughness  can  give  even  greater  preference  to  the  propagation  of  Sn  and  Lg  phase 
velocities  in  the  crust. 

It  may  be  remarked  that  it  is  not  neccessary  for  the  Moho  to  be  rough  at  all  places; 
patches  of  roughness  would  be  sufficient.  Reflection  profiling  of  the  continental  Moho  (e.g., 
Braile  and  Chang,  1986)  suggest  this  could  be  the  case.  The  semblance  analysis  (Figure  2) 
shows  bursts  of  coherent  energy  within  the  P  coda;  these  bursts  could  be  associated  with 
specific  areas  of  roughness. 

The  Role  of  Near-Surface  Intracrustal  Layers 

Can  rough  intracrustal  interfaces  be  partially  responsible  for  peak  Sn  and  Lg  phase  velocities 
observed  in  the  F-K  analysis  of  P  coda?  Figure  4b  gives  an  example  of  a  simplified  crustal 
geometry  which  can  contribute  post-critical  shear  energy  to  the  lower  crust.  In  this  case, 
an  upward  propagating  P  wave,  travelling  in  a  constant  velocity  lower  crust,  scatters  as 
it  encounters  a  highly  irregular  soil-basement  interface.  Incident  energy  is  converted  into 
post-critical  shear  energy  which  is  then  favored  for  lateral  propagation  in  the  crust.  Looking 
at  the  statistical  reflection  coefficients  for  an  irregular  soil-basement  interface,  it  is  apparent 
that  such  near  surface  scattering  may  actually  give  a  far  stronger  contribution  to  post-critical 
shear  energy  than  the  Moho.  Figure  7  shows  reflection  coefficients  for  a  P  wave  incident,  from 
the  lower  crust,  on  a  soil-basement  interface  with  an  rms  slope  of  30°.  It  is  clear  that  this 
interface  can  scatter  the  same  order  of  energy  as  the  Moho  into  phases  which  travel  with  Sn 
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and  Lg  phase  velocities,  and  this  conversion  can  now  occur  not  only  at  grazing  angles,  but 
also  at  most  other  incident  angles,  generating  additional  post-critical  energy  in  the  lower 
crust.  While  the  soil-basement  reflection  case  is  shown  for  illustration  here,  we  generally 
expect  strongly  irregular  interfaces  present  in  the  near  surface  intracrustal  region  (including 
the  free  surface)  to  strongly  influence  the  production  of  scattered  energy,  especially  at  lower 
phase  velocities. 

Just  as  the  scattering  from  the  bottom  of  a  highly  irregular  near  surface  interface  can 
create  large  amplitude,  post-critical  energy  which  remains  trapped  in  the  crust,  the  scattering 
of  P  wave  energy  incident  from  the  top  of  a  near  surface  interface  can  generate  large  amplitude 
energy  which  is  favored  for  lateral  propagation  over  short  distances  in  the  shallower  crustal 
layers,  as  shown  in  Figure  4c.  Figure  7  shows  the  reflection  coefficients  presented  for  the  soil- 
basement  interface,  given  a  P  wave  incident  from  the  soil  layer.  These  reflection  coefficients 
show  that  large  amounts  of  scattered  P  and  converted  SV  energy  is  scattered  to  post-critical 
angles.  This  energy  would  likely  be  observable  only  at  nearby  stations,  since  waves  tend  to 
be  highly  attenuated  in  these  shallower  layers. 

Additional  Remarks  on  Interface  Scattering 

The  above  discussion  assumes  2-D  interface  structure.  The  results  for  the  less  irregular 
Moho  (0  =  10°)  may  be  compared  with  the  3-D  single-scattering  perturbation  results  of 
Prange  (1989).  There  is  a  strong  similarity  in  the  scattering  results  over  interfaces  with 
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similar  Gaussian  statistics.  Post-critical  peaks  and  nulls  match  well  in  both  scattering  angle 
and  relative  amplitudes,  strengthening  the  assertion  that  2-D  scattering  can  give  us  insight 
into  the  mechanisms  of  3-D  scattering.  It  should  be  noted  that  the  reflection  coefficients 
presented  by  Prange  (1989)  do  not  include  specular  contributions. 

Discussion  and  Conclusions 

An  analysis  of  events  recorded  at  the  NORESS,  FINESA,  ARCESS,  and  NYNEX  arrays 
using  F-K  spectra  showed  that  backazimuths  of  P  coda  energy  cluster  around  the  P  backaz- 
imuth  with  a  spread  of  about  -1-/-15°.  The  F-K  spectra  contained  tight  peaks  which  do  not 
spread  over  large  azimuths.  This  contrasts  to  the  behavior  of  S  coda,  and  indicates  that  a 
model  of  scattering  from  a  distributed  homogeneity  is  not  adequate,  at  least  for  the  coherent 
portion  of  the  coda.  The  velocities  measured  fall  into  three  main  groups:  a  high  velocity 
group  (>  7  km/s),  a  Pg  velocity  (6  km/s)  group,  and  a  low  velocity  (<  5  km/s)  group. 
The  high  velocity  group  fits  the  multiple  crustal  reflection  model,  and  is  also  favored  by  P-P 
scattering  from  a  rough  interface.  The  other  velocity  groups  may  be  explained  in  term  of 
scattered,  nonspecular  reflections  from  a  rough  interface.  Scattering  from  an  irregular  Moho 
and  irregular  near  surface  intracrustal  interfaces  are  both  shown  to  be  consistent  with  these 
observations.  It  should  be  noted  that  the  observations  reported  here  were  all  taken  on  stable 
continental  crust,  and  all  conclusions  apply  to  this  situation.  The  more  irregular  the  Moho 
discontinuity  or  the  shallow  intracrustal  interface,  the  more  easily  energy  is  converted  to  the 
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phase  velocities  observed  at  the  arrays.  In  addition,  a  near  surface  interface  in  the  shallow 
crust  can  also  efficiently  generate  scattered  energy  which  is  favored  for  lateral  propagation 
in  the  upper  crustal  layers.  The  effects  of  large  attenuation  in  shallower  crustal  layers  may 
decrease  the  importance  of  these  waves.  This  preliminary  study  indicates  that  a  composite 
model  for  scattering  within  the  crust  must  consider  the  possibility  of  scattering  from  interface 
irregularities. 
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SYN.  RECORD  SECT,,  SMOOTH  INTERFACES 
(TOKSOZ  ET  AL.,  1990) 
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Figure  1;  A  synthetic  seismogram  calculated  with  reflectivity  over  a  simply  layered  crust 
following  a  velocity  model  for  NORESS,  which  is  shown  in  the  text.  The  various  phases 
are  labeled. 
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Figure  2:  The  frequency-wavenumber  spectrum  for  the  P  coda  window  shown.  The  position, 
phase  velocity,  and  azimuth  of  the  primary  peaks  are  given.  Also  shown  is  the  semblance 
analysis  for  this  time  window. 


Figure  3:  Histograms  giving  (a)  the  distribution  of  phase  velocities  and  the  (b)  azimuthal 
deviation  of  events  as  measured  by  hand  from  the  F-K  plots  for  the  data  analyzed  at 
NORESS,  FINESA,  and  ARCESS. 


65 


i 


a) 


c) 


Figure  4;  A  simplified  diagram  showing  (a)  an  irregular  Moho  and  (b)  a  near  surface,  soil- 
basement  interface  generating  energy  which  then  travels  with  Sn-Lg  phase  velocities  in 
P  coda.  Also  shown  is  (c)  a  shallow  near  surface  layer  contributing  to  P  and  S  energy 
which  becomes  trapped  within  a  shallow  surface  layer,  potentially  complicating  Rg  energy 
recorded  at  a  nearby  seismic  array.  The  grayscale  regions  represent  the  corresponding 
mean  reflection  coeflBcient  for  the  given  boundary. 
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Figure  5:  The  numerical  reflection  coefficient  calculated  over  an  irregular  2-D  Moho  discon¬ 
tinuity  with  an  rms  slope  of  10°.  The  total  mean  reflection  coefficient  over  a  randomly 
irregular  interface  with  Gaussian  properties  is  given.  The  S-S  critical  angle  is  approxi¬ 
mately  the  same  as  the  P-P  critical  angle  {9c  =  56°). 


67 


y 


MOHO  REFLECTION  COEFFICIENT 


•90  -70  -50  -30  -10  1  0  30  50  70  90 


scattering  angle  (deg) 


Figure  6:  Similar  to  Figure  5,  except  for  an  interface  with  a  30°. 
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Figure  7:  Numerical  reflection  coefficient  calculated  over  an  irregular  2-D  soil-basement 
interface  with  an  rms  slope  of  30°.  The  total  mean  reflection  coefficient  over  a  randomly 
irregular  interface  with  Gaussian  properties  is  given. 
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